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Abstract 

The unified scaling law for earthquakes, proposed by Bak, Christensen, Danon and 
Scanlon, is shown to hold worldwide, as well as for areas as diverse as Japan, New 
Zealand, Spain or New Madrid. The scaling functions that account for the rescaled 
recurrence-time probability densities show a power- law behavior for long times, with 
a universal exponent about (minus) 2.2. Another decreasing power law governs short 
times, but with an exponent that may change from one area to another. This is in 
contrast with a spatially independent, time-homogenized version of Bak et a/.'s 
procedure, which seems to present a universal scaling behavior. 
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1 Introduction 



Everybody will agree that earthquakes are a complex phenomenon. Indeed, 
one single event can generate a bunch of research papers (several of them 
published even in Nature or Science) , each telling a part of the story of the 
quake. Essentially, these articles would argue that if the tectonic forces were 
like this and that and after the effect of diverse factors and some calculations, 
we would end with precisely that particular earthquake. This research provides 
important information for specific mechanisms triggering earthquakes, "one 
explanation for each earthquake" . 

However, another point of view is possible. As an alternative to the described 
reductionism, Bak was claiming the necessity of a general theory encompass- 
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ing all carthqTiakes [1], from the very small (imperceptible by humans) to the 
largest, catastrophic ones ("killing hundreds of thousands of people"), irre- 
spective of their location at the boundaries or interior of the plates, depth, 
and other tectonic details. A key clue signaling the unity of the phenomenon 
is the existence of the Gutenberg- Richter law [2,3], which states that (for any 
region and for a large enough period of time) the number of earthquakes de- 
creases exponentially as a function of their magnitude. If the seismic processes 
did not form a whole, how could it be that all that variety of events conspire 
together to align onto such a simple curve? 

The first step towards a theory of earthquakes should consist on identifying 
what kind of dynamical process we are dealing with [4]: "is it periodic? Is 
it chaotic? Is it random in space and time?" From our present knowledge, 
the best candidate is self-organized criticahty (SOC) [1,5,6,7]. The analogies 
between earthquakes and SOC systems are clear [8,9,10]: the Earth crust ac- 
cumulates energy (supplied by slow convective motion in the mantle) in the 
form of elastic deformation at a very slow rate. At some point the stress can- 
not be sustained an a rupture initiates, propagating very fast through a fault 
by means of a domino effect, giving rise to an earthquake. So, we have the 
basic ingredients for SOC, i.e., a long-term balance between slow driving and 
fast avalanches in a spatially extended system consisting on many interacting 
parts. The Gutenberg-Richter law is again crucial in this picture, as it im- 
plies that there is not a characteristic scale for the energy dissipated during 
an earthquake (this is so because the energy increases exponentially with the 
magnitude [about a factor 30 in the energy for each unit in the magnitude], 
and therefore, the exponential frequency-magnitude relation transforms into 
a power-law distribution of dissipated energies [3], which is the indication of 
scale invariancc, i.e., criticahty). Therefore, one may talk about the crust as 
being at a critical state, but in contrast to equilibrium critical points, this 
state has to arise spontaneously, as an attractor of the dynamics. 

In addition to the Gutenberg-Richter law, there are other indicators of critical- 
ity or scale invariance in seismicity, as the fractal distribution of hypocenters 
or epicenters [2,3,11], and the Omori law, which tells us that the decay of the 
rate of seismic activity after a large event does not present any characteristic 
time [12]. Further, the structure of faults and tectonic plates is also fractal 
[3,13,14]. For these reasons, although the SOC paradigm is represented by a 
sandpile [15], earthquakes may be considered as the clearest illustration of 
how a real SOC system would look like [1]. 

Despite the initial opposition of the very conservative geophysics community 
to the idea of earthquakes as a SOC phenomenon [1], nowadays SOC is very 
seriously considered by many professional seismologists. It seems that scientific 
evolution in the solid-Earth sciences takes place mostly after great painful 
controversies [16]; and earthquakes are a field where the debate is open [17]. 
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A great hallmark in the carthquake-as-a-SOC-phenomcnon development was 
the release of the Olami, Feder, and Christensen (OFC) model [18] (for a sum- 
mary of the rich dynamical behavior of its variations see Ref. [19]). However, it 
is a subsequent surprising proposal by Ito what calls our interest here [20] : the 
well-studied Bak-Sneppen (BS) model [21], introduced to account for biologi- 
cal evolution, reproduced some other properties of earthquakes. Ito measured 
for a California catalog the same quantities used to characterize the BS model 
[22], in particular, the first- return-time distribution: this is the probability 
that the activity returns (for the first time) at a given spatial location after a 
certain time. As the locations of earthquake occurrence are continuous (in con- 
trast to the BS model), Ito divided the area covered by the catalog into small 
regions of 1° latitude xl° longitude and measured the return times to these 
regions. (Another important difference between the model and real earthquake 
occurrence is that the former is spatially homogeneous, whereas earthquake 
epicenters draw a fractal over the Earth surface.) The results seemed compat- 
ible with a power-law distribution (as in the BS model), but clearly, a more 
in-depth investigation was needed. 

It was Bak, together with Christensen, Danon and Scanlon, who re-opened the 
problem [23,24]. In essence, they used Ito's procedure with the addition of a 
lower bound Mc for the magnitude, in such a way that events with magnitude 
M below Mc were disregarded. This is necessary to avoid spatial and time vari- 
ations in the completeness of the catalogs and to ensure that no events (or not 
many) with M > Mc are missing. A power-law first-return-time distribution 
(followed by a faster decay) was indeed found, but the exponent was different 
from Ito's one (the faster decay was later identified as another power law [25]). 
However, Bak et al. also introduced a crucial element, which was the study of 
the distribution under the variation of the two parameters of the procedure: 
the lower bound Mc and the size L of the small regions. Remarkably, a scaling 
analysis showed that different distributions corresponding to different values of 
L and Mc collapsed onto a single curve under rescaling of the axes by a factor 
XQb^c ji^df ^ where the numerator comes from the Gutenberg- Richter relation 
and the denominator accounts for the fractal distribution of epicenters. In this 
way, it is appropriate to talk about a unified scaling law for earthquakes. 

The law is a direct consequence of Bak's philosophy applied to earthquakes, 
which can be summarized by: 1) Don't care about the tectonic environment. 
(The small regions in which California is divided are independent of it, in con- 
trast with traditional studies.) 2) Don't care about aftershocks, foreshocks, or 
mainshocks (none of these events are removed, all are equally treated, again 
at variance with usual approaches). After all, there is nothing in the seis- 
mograms that differentiates these events; it is only from their relation with 
the other events that these categories can be established, and not without 
unambiguities. 3) Don't care about temporal heterogeneity. 
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We are going to explore in detail Bak et a/.'s unified scaling law; first, intro- 
ducing some variations to their procedure, and next, turning to their method 
and extending it everywhere. Needless to mention, the original flow of ideas 
is the opposite as the one presented here; nevertheless, we believe it is more 
direct to start with our perspective first (at least for us). 



2 Universal scaling law for local distributions of recurrence times 

Let us consider spatial regions of arbitrary shape, which can be of small size 
(as in Ito's paper [20]) but also large. In contrast to Ito and Bak et al. [23,24], 
we concentrate only on one of these regions, where (in the same way as Bak 
et al.) earthquakes with magnitude M above a lower bound Mc are selected. 
Then, if denotes the time of occurrence of the i—th earthquake in the spatial 
and magnitude windows considered, we calculate the recurrence time between 
events i and i — 1 as r, = ij— This time is the same first-return time defined 
by Ito, or the waiting time in the language of Bak et al., or interoccurrence or 
inter-event time in other papers. 

The probability density of the recurrence times can easily be obtained; how- 
ever, in order to pay attention to all the time scales involved in the process it 
will be convenient to look at time in logarithmic scale. The first data set to 
start these measurements is a global earthquake catalog, as the rate of seismic 
occurrence there (defined as number of earthquakes per unit time) is fairly 
constant, and therefore stationary. For instance, for the NEIC worldwide cat- 
alog [26], covering a period from 1973 to 2002 (included), and for several very 
large regions we get the results displayed in Fig. 1 (top curve), after rescaling 
the axes by the seismic rate. All the distributions lie onto a single curve, so. 



where, for a region of spatial coordinates xy, D^yi^r) is the probability density 
that the recurrence time is around a value r and Rxy is the mean rate of seismic 
occurrence (or activity). It is implicit that both Dxy{T) and R^y depend as 
well on Mc and the size of the region; to be concrete, if the region is kept fixed 
the rate depends exponentially on the magnitude, following the Gutenberg- 
Richter relation: R^y = N^ylO'''^'', with N^y the (hypothetical) number of 
events per unit time in the region with magnitude above 0. 

The scaling function / can be represented by 



Dxyi^') Rxyf{RxyT\ 



(1) 




(2) 
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Fig. 1. Single-region recurrence-time probability densities, D^y, after rescaling with 
the rate Rxy The five data sets correspond, from top to bottom, to: 1, the NEIC 
worldwide catalog for regions with L > 180°, 1973 to 2002; 2, NEIC with L < 90° 
(same period of time); 3, Southern California, 1984 to 2001, 1988 to 1991, and 
1995 to 1998; 4, Northern California, 1998 to 2002, 5, Japan, 1995 to 1998, and 
New Zealand, 1996 to 2001; 6 (bottom), Spain, 1993 to 1997, New Madrid, 1975 to 
2002, and Great Britain, 1991 to 2001. The distribution sets have been multiplied 
by 10°, 10-2, 10"^, 10-6, 10-^, and 10-^°, for clarity sake. A total of 82 different 
distributions are shown, with the size of the regions from 0.16° to the whole world, 
and Mc from 1.5 to 7.5. Recurrence times go from 2 min to about 1.5 years; values 
of r < 2 min are not shown. The lines correspond to f{9). 

which has indeed a very general shape. If 7 and S are positive, the former 
controls the shape for small 6 and S the shape at large 6; the situation is 
reversed if both parameters are negative; a is a scale parameter and C a 
normalization correction. 

If regions of smaller size are considered, the rate turns nonstationary, giving 
rise to heterogeneities in time. This is due to large earthquakes, which pro- 
voke a kind of "avalanches of earthquakes", i.e., the aftershock sequences. In 
order to compare with the worldwide case, we concentrate in space-time win- 
dows in which the rate keeps stationary, in other words, we stay away of time 
periods which include very prominent aftershock sequences (by now), in oppo- 
sition to Bak et al. (The simplest way to recognize stationarity is by a linear 
increase with time of the accumulated number of earthquakes in a region.) 
Using the NEIC data and several regional and local catalogs (Southern Cali- 
fornia, Japan, Spain, Great Britain [26], New Zealand, New Madrid [27], and 
also Northern California [28]), we find the results displayed in Fig. 1, taking 
regions of L degrees in longitude and L degrees in latitude. The behavior of 
the distributions is identical to the previous case, collapsing under rescaling 
onto the same universal curve /. (The deviations at short times are due to a 
nonstationary rate at this time scale, provoked by small aftershock sequences.) 
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A massive least-square fit using many regions from tlie NEIC catalog and 
several values of and L yields [26] 7 = 0.67 ±0.05, 5 = 1.05 ±0.05 and a = 
1.64 ±0.15; so, 5 can be considered to be one and we have a decreasing power 
law with exponent 1 — 7 about 0.3 accelerated by an exponential term at large 
times. This type of distribution indicates that earthquakes cluster in time, not 
only for sequences of aftershocks (as it is well known) but even for "background 
seismicity". The counterintuitive consequences of this phenomenon for the 
time evolution of seismic hazard are analyzed in Ref. [27]. 

Finally, the scaling law we propose is valid beyond the stationary limit, re- 
placing the mean rate R^y by its instantaneous value, r^ylt). In this way, the 
probability density ^ of = r,j.y(t)T satisfies ^ = f{rxy{t)T), for aftershock 
sequences where the rate decays following the Omori law, rxy{t) = A^y/t^, 
with the origin of time at the mainshock. Remarkably, the scaling function / 
turns out to be the same as in the stationary case [26]. 

As all the data analyzed are well fit by the same function /, with the same 
values of the parameters, we may talk about a universal scaling law. Neverthe- 
less, this law cannot be designated as unified, since the scaling only includes 
the Gutenberg- Richtcr law, but not the fractal dimension of the epicenters, in 
contrast to Bak et al.^s law. And in terms of its definition in space, D^y may 
be referred to as a local distribution (although the size of the xy— region can 
grow to reach the total area covered by the catalog) . 



3 Unified scaling law beyond Southern California 

Now it is time to pay attention to the original approach of Bak et al. We 
only have to consider all the (non-overlapping) regions of size L necessary to 
cover completely a much larger area (Southern California, Japan, the whole 
world, or New Madrid) and expand the time windows without bothering about 
the nonstationarities of the seismic rate. We measure for each a;|/— region the 
recurrence times r/^^^ in the same way as before, with the difference that all 
these series of recurrence times are counted into one single probability density, 
T>{t), which we may call global (in contrast with our local version). It is found 
that V^t) scales, under the change of L and the lower bound Mc, as: 

V(r) ^nF(nr), (3) 

with TZ the mean value of the local mean rate R^y, calculated over all the 
xy— regions of size L with seismic activity, i.e., TZ = J2xy Rxy/n, where n is 
the number of such regions. From the equation for R^y and the scaling of n 
with L, n = {e/iyf, we get, 7^ = N{L/eyf 10"^^% with = T^xy N^y and £ 
a rough measure of the linear size of the total area under study (in degrees) . 
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Fig. 2. Recurrence-time probability densities, calculated following Bak et aZ.'s pro- 
cedure, after rescaling by TZ. 84 distributions are shown, with L ranging from 0.039° 
to 45°, and 1.5 < Mc < 6. The curves are shifted by factors 10°, 10"^ 10"^ and lO'^, 
and correspond to: 1 (top), Southern California, 1984 to 2001; 2, Northern Califor- 
nia, 1985 to 2003; 3, Southern Cahfornia, 1988 to 1991 (stationary rate), NEIC, 1973 
to 2002, Japan, 1995 to 1998, and Spain, 1993 to 1997; 4 (bottom). New Zealand, 
1996 to 2001, and New Madrid, 1975 to 2002. Short or intermediate times are fit 
by: 1 (top), 0.12/6'°-9^; 2, 0.15/6lO-9, 3 and 4 (bottom), 0.05/6'°-95 and 0.5/6'°•^ with 
9 = TZt. In all cases the long-time tail is fit by 0.25/0^'^. The times in the horizontal 
axis span from 2 min to about 20 years. Recurrence times smaller than 2 minutes 
are not shown, except for Japan and the NEIC catalog where the minimum times 
are 4 min and 10 min, respectively. 

In addition to Southern California [23,24], the scaling relation for T) is valid for 
the catalogs studied in the previous section, see Fig. 2. However, the scaling 
function F does not seem to be universal; the clear decreasing two-power-law 
behavior of Southern California [25], with exponent for small r about 0.95 and 
2.2 for large r, changes slightly for Northern California (exponents 0.9 and 2.2), 
but becomes more complicated for the other catalogs, with the appearance of 
an intermediate bump. Figure 2 shows how this bump can be approximated 
by another decreasing power law, with exponent 0.5, roughly. Further, for 
New Madrid there is no trace of the 0.95 exponent, and for New Zealand the 
behavior is not clear. The only exponent that seems to be universal is the 
one for large times, 2.2 in all cases. The unified scaling law could also hold 
for Great Britain, but the few data considered there, only about 500 events, 
makes the statistics too poor for small and intermediate times. 

Let us mention that the probability densities of Bak et al. and ours, in addition 
to quite different shapes, have different meanings too. In our case, we provide 
the probability of return of an earthquake for a given xy— region of size L, 
with the only information required being the mean seismic rate there, R^y 
On the other hand, Bak et al.'s distribution gives the return probability if one 
does not know in which region of size L of a much larger area (like Southern 
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California) one is (or you know the region but you don't have knowledge about 
the rate Rxy), and it is necessary to know the average rate TZ for regions of 
size L. The relation between these two approaches are studied in Ref. [25]. 

As a final conclusion, we see that scaling is an intrinsic characteristic of seis- 
micity. Both the unified scaling law and our local approach show that the 
distributions of recurrence times scale with some average value of the rate 
of seismic activity. This is a clear consequence of the scale-invariant struc- 
ture of seismic occurrence in time, space and magnitude. Whatever triggers 
earthquakes operates in the same way at all spatial and temporal scales. 
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